\documentclass[a4paper,10pt]{article}
\usepackage{amsmath,amsfonts,amsthm,amssymb,graphicx, url}
\usepackage{listings}
\lstset{language=C}
\lstset{basicstyle=\small}
\newcommand{\bvec}[1]{\boldsymbol{#1}} % vectors in bold instead of with an arrow on top

% Title Page
\title{Assignment series 2 -- Computational Finance}
\author{Florian Speelman \& Jannis Teunissen}


\begin{document}
\maketitle
\section{Stock price model}
We model the evolution of the stock price $S$ as a geometric Brownian motion:
\[
 dS = r S dt + \sigma_0 S dz,
\]
where $r$ is the risk free return, $\sigma_0$ is the volatility and $z$ is a Wiener process.
This means that $dz = \epsilon \sqrt{dt}$, where $\epsilon$ is a random drawing from a standard
normal distribution.

It\^{o}'s Lemma states that a function $G(S,t)$ will changes in time as
\[
 dG = \left[\frac{\partial G}{\partial S}r S + \frac{\partial G}{\partial t}
+\tfrac{1}{2}\frac{\partial^2 G}{\partial^2 S}\sigma_0^2 S^2\right]dt
+ \frac{\partial G}{\partial S}\sigma_0 S dz.
\]
If we take $G(S,t) = \ln(S)$ we thus get
\[
 d\ln(S) = (r - \tfrac{1}{2}\sigma_0^2)dt + \sigma_0 dz.
\]
We assume $r$ and $\sigma_0$ to be constant, so that $\ln(S)$ is a generalized Wiener process.
Now if $S_T$ denotes the stock price
at time $T$, then we know that $\ln(S_T)$ will have a normal distribution with mean
\[
 E[\ln(S_T) - \ln(S_0)] = (r - \tfrac{1}{2}\sigma_0^2)T,
\]
and variance
\[
 \operatorname{Var}[\ln(S_T) - \ln(S_0)] = \sigma_0^2 T.
\]
If $\Phi(\mu,\sigma)$ denotes a normal distribution with mean $\mu$ and standard deviation $\sigma$,
we can express the distribution of $S_T$ as
\[
  S_T \sim \exp\left[\ln(S_0) + \Phi\left((r - \tfrac{1}{2}\sigma_0^2)T,\sigma_0 \sqrt{T}\right)\right],
\]
or the equivalent form
\[
 S_T \sim S_0\exp\left[(r - \tfrac{1}{2}\sigma_0^2)T + \sigma_0\sqrt{T}\Phi(0,1)\right].
\]
With the above formula it is straightforward to sample from $S_T$: draw a random
number $\epsilon$ from a standard normal distribution, then the sample $s_T$ is given by
\begin{equation}
s_T = S_0\exp\left[(r - \tfrac{1}{2}\sigma_0^2)T + \sigma_0\sqrt{T}\cdot\epsilon\right] 
\label{eq:stsample}
\end{equation}

\section{Monte Carlo simulation statistics\label{statintro}}
Consider a random variable $\chi$ with mean $\theta$ and variance $\sigma^2$. If we take the average of $n$ independent
samples $\chi_i$ we have an estimate $\bar{\chi}$ of the mean
\[
 \bar{\chi} = \frac{1}{n}\sum_{i=1}^{n}\chi_i.
\]
Since the $\chi_i$ are independent, the variance of $\bar{\chi}$ is $\sigma^2/n$. If the number of samples $n$ is large,
by the central limit theorem $\bar{\chi}$ will approximately be normally distributed. Then we can estimate the probability of
the sample mean lying within $a$ standard deviations as
\[
 P\left(|\bar{\chi} - \theta| < \frac{a\sigma}{\sqrt{n}}\right) = \operatorname{N_{CDF}}(a) - \operatorname{N_{CDF}}(-a),
\]
where $\operatorname{N_{CDF}}$ is the cumulative standard normal distribution function.
Usually the variance $\sigma^2$ is not known, but for a large number of samples we can use the unbiased
estimator $\hat{\sigma}^2$:
\[
 \hat{\sigma}^2 = \frac{1}{n-1}\sum_{i=1}^{n}(\chi_i - \bar{\chi})^2.
\]
The standard error of $\bar{\chi}$ is given by $\hat{\sigma}/\sqrt{n}$, so that for an $k$ times smaller
standard error we need about $k^2$ times as many samples.

Typically the results of a Monte Carlo simulation will look like figure \ref{fig:monteresult}, where every
point was generated using a different seed for the randum number generator. It should be clear that most
points lie within a few standard errors of the mean of the distribution. See figure \ref{fig:conveuc} for
an illustration of the fact that the standard error decreases like $1/\sqrt{n}$.
\begin{figure}
\begin{center}
 \includegraphics[width = 9cm]{conv-eu-call.png}
\caption{
Every point shown is the result of a Monte Carlo simulation with a different seed and a different number of trials,
for an European call option with parameters as in equation \eqref{eq:standardpars}. The plotted error bars have a length
of two standard errors, and we see that all points lie within a few standard errors of the mean. 
\label{fig:monteresult}}
\end{center}
\end{figure}
\subsection{Antithetic variables}
The variance of the sum of two  random variables  $\chi_i + \chi_{i+1}$ is given by
\[
 \operatorname{Var}(\chi_i + \chi_{i+1}) = \operatorname{Var}(\chi_i) + \operatorname{Var}(\chi_{i+1})
 + 2\operatorname{Cov}(\chi_i, \chi_{i+1}).
\]
If $\chi_i$ and $\chi_{i+1}$ are negatively correlated, we can decrease the variance in our simulation
by taking the average of $n/2$ pairs $\chi_i + \chi_{i+1}$ instead of $n$ independent samples. We try to value
options in our simulation, and each sample $\chi_i$ of the option value is based on $M$ drawings
$\epsilon_1, \dots, \epsilon_M$ from a standard normal distribution:
\[
 \chi_i = f(\epsilon_1, \dots, \epsilon_M).
\]
It can be proven that if $f$ is a monotone function in all of its arguments, then it will be negatively
correlated to
\[
 \chi_{i+1} = f(-\epsilon_1, \dots, -\epsilon_M).
\]
We can not only get a smaller variance (or standard error) by using pairs $\chi_i + \chi_{i+1}$, but we also
have to generate only half the random numbers compared to $n$ independent samples.

However, it will depend on the option how negative the correlation is, and it can even become positive.
If the mean of $\chi_i$ is $\theta$, we can express the covariance as
\[
 \operatorname{Cov}(\chi_i, \chi_{i+1}) = E[(\chi_i - \theta)(\chi_{i+1} - \theta)].
\]
For an European option this can never become positive (loosely speaking, if one goes down, the other
goes up, so the covariance will be negative). But if we for example have a portfolio of two barrier options,
one of type up-and-out and the other of type down-and-out, the covariance can be positive. Then the use
of antithetic variables leads to a greater variance and standard error.

As an example consider an European call option with the following properties
\begin{equation}
 S_0 = 100,\;\;K = 99,\;\;\sigma_0 = 0.2,\;\;r = 0.06,\;\;T = 1,
\label{eq:standardpars}
\end{equation}
where $S_0$ is the current stock price, $K$ is the strike, $\sigma_0$ is the volatility, $r$ the risk free
return and $T$ the time to maturity. Figure \ref{fig:conveuc} shows how the standard error
decreased with the number of trials in a simulation, with and without antithetic variables.
\begin{figure}
\begin{center}
 \includegraphics[width = 9cm]{conv-eu-call-stddev.png}
\caption{Standard error in the value of an European call option, with parameters as in equation \eqref{eq:standardpars}.
Both curves are proportional to $1/\sqrt{n}$, where $n$ is the number of trials,
but the use of antithetic variables led to a reduction of the standard error of about $40\%$.
\label{fig:conveuc}}
\end{center}
\end{figure}
\subsubsection{Variance and covariance of an European call option}
To estimate the variance reduction that can be obtained using antithetic variables, we need to know both the variance
of $\chi_i$ and the covariance of a antithetic pair $\chi_i, \chi_{i+1}$. In this section we will derive expressions
for both these terms that may also be useful for using a European call option as a control variate.

Consider an European call option with time to maturity $T$ and strike price $K$.
The value of
a particular sample $\chi_i$ is $e^{-rT}(s_T - K)^{+}$, where $s_T$ is given by equation \eqref{eq:stsample}:
\begin{equation*}
 \chi_i = e^{-rT}\left(S_0\exp\left[(r - \tfrac{1}{2}\sigma_0^2)T + \sigma_0\sqrt{T}\epsilon\right] -K\right)^{+}.
\label{eq:chi}
\end{equation*}
The variance of $\chi_i$ is given by
\begin{align*}
 \operatorname{Var}(\chi_i) &= E(\chi_i^2) - E(\chi_i)^2\\
 &= e^{-2rT}E\left[((S_T - K)^{+})^2\right] - c^2\\
 &= e^{-2rT}\int_{K}^{\infty}(S^2 -2 K S + K^2)f\left(S,\mu,\sigma\right)dS - c^2.
\end{align*}
The last term is the square of the Black-Scholes value $c$ of this option, 
and $f(S,\mu,\sigma)$ is the probability density function of $S_T$:
\begin{equation}
 f(S, \mu, \sigma) = \left(S \sigma \sqrt{2\pi}\right)^{-1} \exp\left[-\left(\ln S - \mu\right)^2/\left(2\sigma^2\right)\right].
\label{eq:probdens}
\end{equation}
The approriate values for $\mu$ and $\sigma$ are $\ln(S_0) +(r - \tfrac{1}{2}\sigma_0^2)T$ and $\sigma_0 \sqrt{T}$ respectively.
In our previous report we have already calculated integrals like the above (section 1.7), and from the results we got there
it is straightforward to derive
\begin{align*}
 \int_{K}^{\infty}2 K S f\left(S,\mu,\sigma\right)dS
 &= 2 K S_0 e^{rT} \operatorname{N_{CDF}}(d_1),\\
 \int_{K}^{\infty}K^2 f\left(S,\mu,\sigma\right)dS
 &= K^2 \operatorname{N_{CDF}}(d_2),\\
 \int_{K}^{\infty}S^2 f\left(S,\mu,\sigma\right)dS
 &= S_0^2\exp\left[(2r + \sigma_0^2)T\right] \operatorname{N_{CDF}}(d_3),
\end{align*}
where the $d_i$ are given by
\begin{equation*}
 d_1 = \frac{\ln(S_0/K) + (r + \sigma_0^2/2)T}{\sigma_0\sqrt{T}},\;\; d_2 = d_1 - \sigma_0\sqrt{T},\;\; d_3 = d_1 + \sigma_0\sqrt{T}.
\label{d1d2}
\end{equation*}
So the variance is equal to
\begin{align}
 \operatorname{Var}(\chi_i) =& 2 e^{-rT} K S_0 \operatorname{N_{CDF}}(d_1)[1+\operatorname{N_{CDF}}(d_2)]\nonumber\\
 &+ e^{-2rT}K^2 \operatorname{N_{CDF}}(d_2)[1-\operatorname{N_{CDF}}(d_2)]\nonumber\\
 &- S_0^2 \operatorname{N_{CDF}}(d_1)^2 + S_0^2 e^{\sigma_0^2T}\operatorname{N_{CDF}}(d_3).
\label{eq:var}
\end{align}
If we use antithetic variables, the covariance between $\chi_i$ and $\chi_{i+1}$ is given by
\begin{align*}
 \operatorname{Cov}(\chi_i, \chi_{i+1}) &= E(\chi_i \cdot \chi_{i+1}) - E(\chi_i)E(\chi_{i+1})\\
& = e^{-2rT}E\left[(s_T - K)^{+}(\hat{s}_T - K)^{+}\right] - c^2,
\end{align*}
where $\hat{s}_T$ is obtained by using $-\epsilon$ instead of $\epsilon$ in equation \eqref{eq:stsample}.
Now since $\hat{s}_T = 1/s_T\cdot S_0^2 e^{(2r-\sigma_0^2)T}$, for both of the terms with the $^{+}$
to be positive we need to have $s_T > K$ and $s_T < \beta / K$, where $\beta = S_0^2 e^{(2r-\sigma_0^2)T}$.
Thus the covariance can be expressed as
\begin{align}
 \operatorname{Cov}(\chi_i, \chi_{i+1}) = e^{-2rT}\int_{K}^{\beta/K}\left[\beta -K(S + \beta/S) + K^2\right]f\left(S,\mu,\sigma\right)dS
 - c^2,
\label{eq:covformula}
\end{align}
where the integral must be taken zero when $\beta/K < K$. 
The above formula can be worked out like we did for the variance, but we won't do that here.
If we had some more time it would be interesting to compare the above formulas to experiment, but
for some discussion see figures \ref{fig:sevsstrike} -- \ref{fig:sevssigma}.
\begin{figure}
\begin{center}
 \includegraphics[width = 9cm]{conv-strike.png}
\caption{Standard error in the value of an European call option, with parameters as in equation \eqref{eq:standardpars}
and one million trials. Clearly antithetic variables give less of an improvement for $K \gtrsim S_0$. If we look at equation
\eqref{eq:covformula} this is not so strange, for `large' $K$ both the integral and $c$ (the Black-Scholes
value of the option) will be small compared to the variance (equation \eqref{eq:var}). However, for a low strike price
the use of antithetic variables leads to a significant reduction of the standard error.
\label{fig:sevsstrike}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
 \includegraphics[width = 9cm]{conv-sigma.png}
\caption{Standard error in the value of an European call option, with parameters as in equation \eqref{eq:standardpars}
and two million trials. With and without antithetic variables the standard error increases with $\sigma_0$, but the difference
in standard error between the two methods remains almost constant.
\label{fig:sevssigma}}
\end{center}
\end{figure}


\subsection{Control variates\label{controlvars}}
In our Monte Carlo simulation we are trying to estimate the mean $\theta$ of $\chi$.
Suppose there exists some other random variable $Y$ with known mean $\theta_Y$ that is correlated to $\chi$,
then
\begin{equation*}
  \chi_c = \chi + \alpha(Y - \theta_Y)
\label{eq:ctrlvar}
\end{equation*}
is an unbiased estimator of $\theta$ for any value of $\alpha$. The variance of $\chi_c$
is given by
\[
 \operatorname{Var}(\chi_c) = 
\operatorname{Var}(\chi) +  \alpha^2 \operatorname{Var}(Y) + 2\alpha\operatorname{Cov}(\chi, Y).
\]
This is minimized when $\alpha = -\operatorname{Cov}(\chi, Y)/\operatorname{Var}(Y)$, and if we divide
by $\operatorname{Var}(\chi)$ we get
\[
 \frac{\operatorname{Var}(\chi_c)}{\operatorname{Var}(\chi)} = 
1 - \operatorname{Corr}(\chi,Y)^2.
\]
So the variance will be greatly reduced if the correlation is close to one. Usually both $\operatorname{Var}(Y)$
and $\operatorname{Cov}(\chi, Y)$ are not known, but these can be estimated from the simulation data:
\[
 \operatorname{Var}(Y) \approx \frac{1}{n-1}\sum_{i=1}^{n}(Y_i - \bar{Y})^2,
\]
\[
 \operatorname{Cov}(\chi, Y) \approx \frac{1}{n-1}\sum_{i=1}^{n}(\chi_i - \bar{\chi})(Y_i - \bar{Y}).
\]

Whether the use of a control variate pays off depends on the correlation, but also on the cost of simulating the
control variate. For example, if an European option is used as a control variate for an American option,
then the simulation of the European option is quite cheap, since there is no path-dependence.

\section{Value of a geometric average Asian option\label{asiang}}
A geometric average Asian call option with strikeprice $K$ and time to maturity $T$ has an expected value $c$
\[
 c = e^{-rT} E[(\tilde{S} - K)^{+}],
\]
where the tilde denotes the geometric average: $\tilde{S} = \left(\prod_{i=1}^{N}S_i\right)^{1/N}$.
Here $S_i$ is the stockprice at time $t_i = \tfrac{i}{N} T$. Given some $S_i$, $S_{i+1}$ is distributed as
\begin{equation}
\ln\left(\frac{S_{i+1}}{S_i}\right) \sim \Phi\left[(r - \sigma^2/2)\tfrac{T}{N}, \sigma\sqrt{T/N}\right],
\label{eq:siplusone}
\end{equation}
where $\Phi(\mu,\sigma)$ is a normal distribution with mean $\mu$ and standard deviation $\sigma$. Let
us rewrite $\tilde{S}$ in the following way
\[
 \tilde{S} = \left(\frac{S_N}{S_{N-1}}\right)^{1/N}
	      \left(\frac{S_{N-1}}{S_{N-2}}\right)^{2/N}
	      \left(\frac{S_{N-2}}{S_{N-3}}\right)^{3/N}
	      \dots
	      \left(\frac{S_{1}}{S_{0}}\right)^{N/N} S_0.
\]
If we now use equation \eqref{eq:siplusone} we see that $\ln(\tilde{S})$ is distributed as
\[
 \ln(\tilde{S}) \sim \ln(S_0) + \sum_{i=1}^{N}\frac{i}{N}\cdot\Phi\left[(r - \sigma^2/2)\tfrac{T}{N}, \sigma\sqrt{T/N}\right].
\]
The sum of two indepedent normal random variables $\alpha\chi_1$ and $\beta\chi_2$ is again a normal random variable,
with mean $\alpha E(\chi_1) + \beta E(\chi_2)$ and variance
$\alpha^2 \operatorname{Var}(\chi_1) + \beta^2 \operatorname{Var}(\chi_1)$. So $\ln(\tilde{S})$ is a normal random
variable with mean $\tilde{\mu}$
\[
\tilde{\mu} = \ln(S_0) + \frac{N+1}{2N}\cdot(r - \sigma^2/2)T,
\]
and variance $\tilde{\sigma}^2$
\[
 \tilde{\sigma}^2 = \frac{(N+1)(2N+1)}{6N^2}\sigma^2 T.
\]
This is very similar to an European style option, only the mean and variance of the random variable are different.
Therefore, the Black-Scholes value of this Asian option is given by
\footnote{See section 1.7 of our previous report for the derivation of the Black-Scholes analytical value
of an European style option}
\begin{align}
 c = 
\exp\left[\tilde{\mu} + \tilde{\sigma}^2/2 - rT\right] &\cdot \operatorname{N_{CDF}}\left(\frac{-\ln K +\tilde{\mu} + \tilde{\sigma}^2}{\tilde{\sigma}}\right)\nonumber\\
&- K e^{-r T} \cdot \operatorname{N_{CDF}}\left(\frac{-\ln K +\tilde{\mu}}{\tilde{\sigma}}\right),
\label{eq:cvalue}
\end{align}
Introducing two new variables
\begin{align*}
 \tilde{d}_1 &= \sqrt{\frac{6N^2}{(N+1)(2N+1)}} \cdot
  \frac{\ln(S_0/K) + \frac{N+1}{2N}rT + \frac{(N+1)(N+2)}{12N^2}\sigma^2 T}{\sigma\sqrt{T}},\\
 \tilde{d}_2 &= \tilde{d}_1 - \sqrt{\frac{(N+1)(2N+1)}{6N^2}}\cdot\sigma\sqrt{T},
\end{align*}
we can rewrite equation \eqref{eq:cvalue}:
\[
 c = S_0 \operatorname{N_{CDF}}(\tilde{d}_1) \exp\left[\frac{1-N}{2N}rT + \frac{1-N^2}{12N^2}\sigma^2T\right]
   - K e^{-rT} \operatorname{N_{CDF}}(\tilde{d}_2).
\]
Note that for $N = 1$ the above expressions simplify to a well known result: the Black-Scholes value of
an European style option. This is to be expected, since an American option is equivalent to an Asian option that
is evaluated only at $t = T$ (so for which $N = 1$).

The above derivation was done for a call option, but the value of a put option can be derived in a very similar way.
\subsection{Simulating an arithmetic Asian option}
An arithmetic asian option is very similar to a geometric one, except for that it uses the arithmetric average. It
is hard to derive an analytical expression for the value of such an option.

Using the results from section \ref{controlvars} and \ref{asiang} we wrote a program that could perform
Monte Carlo simulations using control variates.
For a better comparison we consider a sample with a control variate as two trials, since basically two trials are done:
one for the random variable itself and one for the control variate.
See figures \ref{fig:asianstrike}--\ref{fig:asianstdvstrials} for results and discussion.

\begin{figure}
\begin{center}
 \includegraphics[width = 9cm]{asian-strike.png}
\caption{
Standard error for an arithmetic Asian option with parameters as in equation \eqref{eq:standardpars}, $N = 365$ and
$2\cdot10^5$ trials. We see the standard error decrease for large $K$ because the option value goes to zero.
The results in this figure are from simulations without a control variate, compare with figure \ref{fig:asianctrstrike}.
\label{fig:asianstrike}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
 \includegraphics[width = 9cm]{asian-strike-cv.png}
\caption{
Standard error for the same option as in figure \ref{fig:asianstrike} but now we have used a geometric Asian option as control variate.
The standard error is about 25 times smaller using the control variate, and behaves in a similar way as before
for increasing $K$. This is quite a significant improvement, withouth a control variate about 600 times as many trials would
be necessary to achieve this.
\label{fig:asianctrstrike}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
 \includegraphics[width = 9cm]{asian-points-cv.png}
\caption{
Standard error for the same option as in figure \ref{fig:asianstrike} for different $N$. Without a control
variate the standard error is a lot larger, but the behaviour with increasing $N$ is very similar: in both
cases the standard error seems to be largely indepedent of $N$.
\label{fig:asianstdvsn}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
 \includegraphics[width = 9cm]{asian-trials-cv.png}
\caption{
Standard error for the same option as in figure \ref{fig:asianstrike} for a different number of trials.
With or without a control variate we see the same behaviour: the standard error decreases as
$1/\sqrt{n}$ where $n$ is the number of trials. This is in agreement with section \ref{statintro},
since we still look at the average of $n/2$ independent samples.
\label{fig:asianstdvstrials}}
\end{center}
\end{figure}

\section{Estimating the hedge parameter}
The hedge parameter $\delta$, or Delta, measures the sensitivity of
an option to the value of its underlying asset. In our case, it is the derivative
of the value $V$ of the option to the stock price $S$:
\[
\delta = \frac{\partial V}{\partial S}\text{.}
\]
We use two different ways of estimating this Greek, bump and revalue
and the likelihood ratio method.

\subsection{Bump and revalue}
The first way to compute $\delta$ we use is the Euler scheme,
in this context also called bump and revalue. The
Monte Carlo simulation is executed twice, once with the original
$S_0$ and once with one increased slightly, by an amount $\epsilon$.
\[
\delta = \frac{V(S_0 + \epsilon) - V(S_0)}{\epsilon}
\]
This gives a variance of
\[
\operatorname{Var}(\delta) = \frac{1}{\epsilon^2}\left[\operatorname{Var}\left(V(S_0 + \epsilon)\right)
				+\operatorname{Var}\left(V(S_0)\right) - 2\operatorname{Cov}\left(V(S_0 + \epsilon), V(S_0)\right)\right]\text{.}
\]
To control the variance we use the same random seed for the bumped and unbumped estimate of the value,
because the covariance then typically is quite large, reducing $\operatorname{Var}(\delta)$, see table \ref{table1}.
\begin{table}
\begin{tabular}{|c|c|c|c|c|}
\hline
&\multicolumn{2}{|c|}{Different seed} & \multicolumn{2}{|c|}{Same seed}  \\ \hline
& $\epsilon = 0.01$ & $\epsilon = 0.5$ & $\epsilon = 0.01$ & $\epsilon = 0.5$\\ \hline
European & $7.411$ & $0.149$ & $1.818 \cdot 10^{-3}$ & $1.817\cdot 10^{-3}$\\ \hline
Digital & $0.218$ & $4.44 \cdot 10^{-3}$ & $4.213 \cdot 10^{-3}$ & $6.110\cdot 10^{-4}$\\ \hline
\end{tabular}
\caption{The standard error in the value of $\delta$, for an option with parameters as in equation \eqref{eq:standardpars}
and $2\cdot10^5$ trials. These standard errors were obtained by comparing 200 full simulations, and it is quite clear
that there is little to no convergence when using different seeds. For these parameters the analytical $\delta \approx 0.6737$.
\label{table1}}
\end{table}

\subsection{Likelihood ratio method}
We also implemented the likelihood ratio method the estimate $\delta$ in our simulations.
In this method we use Monte Carlo on an estimator for the $\delta$.
If we have 
\[
s_T = S_0\exp\left[(r - \tfrac{1}{2}\sigma_0^2)T + \sigma_0\sqrt{T}\cdot Z \right] 
\]
(like in equation (\ref{eq:stsample}) but using $Z$ for the random
variable to avoid confusion with the bump and revalue method)
then
\[
e^{-r T} f(s_T) \frac{Z}{\sigma_0 S_0 \sqrt{T}}
\]
is an unbiased estimator for $\delta$. Here $f(s_T)$ is a function
describing the pay off of the option at maturity. The derivation
of this method is outside of the scope of this report.

For results see table \ref{table2}. This method is easier
to use and analyse than the Euler method, even though the used equations
are harder to derive. 

\begin{table}
\begin{tabular}{|c|c|c|c|c|}
\hline
&\multicolumn{2}{|c|}{European} & \multicolumn{2}{|c|}{Digital}  \\ \hline
Trials & Result & Standard error & Result & Standard error \\ \hline
$10^4$ & $0.668432$ & $0.015592$ & $0.018342$ & $0.000296$ \\ \hline
$10^5$ & $0.681221$	& $0.005126$ & $0.018338$ & $9.414941 \cdot 10^{-5}$ \\ \hline
$10^6$ & $0.672764$ & $0.001602$ & $0.018215$ & $2.965956 \cdot 10^{-5}$ \\ \hline
$10^7$ & $0.674040$ & $0.000508$ & $0.018218$ & $9.389574 \cdot 10^{-6}$ \\ \hline
\end{tabular}
\caption{Results for the likelihood ratio method for both a european
and a digital option. Convergence
goes with the square root of the amount of trials, as is common
with Monte Carlo methods. 
\label{table2}}
\end{table}

\section{Conclusion}
What we have clearly seen is that Monte Carlo simulations are easy to implement, but the standard error converges
quite slow as $1/\sqrt{n}$, where $n$ is the number of trials. It is essential to use variance reduction techniques, and
we did so quite successfully. Both the use of control variates and antithetic variables can lead to very large variance reductions,
depending on the options that are being priced. Unfortunately we did not have time to implement other techniques such as importance sampling,
or the use of pseudo random number generators.

The estimation of the hedge from a Monte Carlo simulation is not very efficient, requiring a whole series of complete simulations.
However, using the bump-and-revalue or likelihood ratio method an accurate estimate is usually possible.
\end{document}
